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NEUMANN-NEUMANN WAVEFORM RELAXATION ALGORITHM 
IN MULTIPLE SUBDOMAINS FOR HYPERBOLIC PROBLEMS IN 

ID AND 2D 


BANKIM C. MANDAL * 

Abstract. We present a Waveform Relaxation (WR) version of the Neumann-Neumann algo¬ 
rithm for the wave equation in space-time. The method is based on a non-overlapping spatial domain 
decomposition, and the iteration involves subdomain solves in space-time with corresponding inter¬ 
face condition, followed by a correction step. Using a Fourier-Laplace transform argument, for a 
particular relaxation parameter, we prove convergence of the algorithm in a finite number of steps 
for finite time intervals. The number of steps depends on the size of the subdomains and the time 
window length on which the algorithm is employed. We illustrate the performance of the algorithm 
with numerical results, followed by a comparison with classical and optimized Schwarz WR methods. 

Key words. Neumann-Neumann, Waveform Relaxation, Wave equation, Domain Decomposi¬ 
tion. 


1. Introduction. We formulate a new variant of Waveform Relaxation (WR) 
methods based on the Neumann-Neumann algorithm to solve hyperbolic problems in 
parallel computer, and present convergence results for the method. The Neumann- 
Neumann algorithm was introduced for solving elliptic problems by Bourgat et al. [1], 
see also [2G] and [28]. The iteration involves solving the subdomain problems using 
Dirichlet interface conditions in the first step, followed by a correction step involving 
Neumann interface conditions. The convergence behavior of the algorithm is now well 
understood for elliptic problems, see for example the book [29]. 

To solve time-dependent problems in parallel, the following three possible classes 
of domain decomposition techniques exist: 

• this approach consists of discretizing the problem uniformly in time with an 
implicit scheme to obtain a sequence of elliptic problems, which are then 
solved by DD methods. For this kind of technique, we refer to [2, 8]. One 
disadvantage of this approach is that, uniform time step across the whole do¬ 
main need to be enforced, which is very restrictive for problems with variable 
coefficients or multiple time scales. Also this method is expensive for parallel 
computation, since one needs to exchange information at each time step of 
the discretization. 

• in this approach the equation is first discretized in space, which is called the 
method of lines, and then one applies a waveform relaxation algorithm to solve 
the large system of ordinary differential equations (ODEs) obtained from the 
space-discretization process. Multigrid dynamic iteration [21, 16] and multi¬ 
splitting algorithms [17] are some particular examples of this approach. 

• in contrast to the two classical techniques above, there exist space-time do¬ 
main decomposition methods, formulated at the continuous level. Here, in¬ 
stead of discretizing in time or in space, one decomposes the original spatial 
domain into smaller subdomains and considers each subproblem as posed in 
both space and time; then the subproblems are solved iteratively communi¬ 
cating information at the interfaces between subdomains. This permits the 
use of different numerical methods in different subdomains. At each itera¬ 
tion, one solves the space-time subproblem over the entire time interval of 
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interest, before communicating interface data across subdomains. Thus one 
saves communication time while computing in parallel computer. For this ap¬ 
proach, see [14, 1.5, 4, 6, 24] for parabolic problems, and [8, 5, 7] for hyperbolic 
problems. 

In this article, we focus on the WR-type algorithm as it allows different discretiza¬ 
tions in different space-time subdomains. WR methods have their origin in the work 
of Picard [2-5] and Lindelof [20] in the late 19th century. Lelarasmee, Ruehli and 
Sangiovanni-Vincentelli [19] rediscovered WR as a parallel method for the solution of 
ODEs. 

In a different viewpoint, the WR-type methods can be seen as an extension of 
DD methods for elliptic PDEs. The systematic extension of the classical Schwarz 
method to time-dependent parabolic problems was started in [14, 15]; later optimized 
SWR methods have been introduced to achieve faster convergence or convergence 
with no overlap, see [6] for parabolic problems, and [8] for hyperbolic problems. 
Recently, we extended the substructuring methods, namely the Dirichlet-Neumann 
and Neumann-Neumann methods, to space-time problems; for parabolic problems see 
[13, 23, 18, 11, 22], and for hyperbolic problems see [12, 11, 22]. We analyzed for the 
heat equation to prove that on finite time intervals, the Dirichlet-Neumann Waveform 
Relaxation (DNWR) and the Neumann-Neumann Waveform Relaxation (NNWR) 
methods converge superlinearly for an optimal choice of the relaxation parameter. 
On the contrary for the wave equation, these methods with a two-subdomains decom¬ 
position converge in a finite no of steps, see [11]. In this paper, we propose the NNWR 
method with many subdomains decomposition for hyperbolic problems and analyze 
the method for the Wave equation in one and two space dimensions. We analyze 
the method in the continuous setting to ensure the understanding of the asymptotic 
behavior of the methods in the case of fine grids. 

We consider the following hyperbolic PDE on a bounded domain D C 0 < t < 
T, d = 1, 2, 3, with a smooth boundary as our guiding example, 


c^{x)Au 

= 

X € il,0 <t < T, 

u{x, 0) 

= Uoix), 

X G 

ut{x,0) 

= Voix), 

X G il, 

u{x,t) 

= 9{x,t), 

X G dfl, 0 <t <T, 


with c{x) being a positive function. 

We introduce in Section 2 the non-overlapping NNWR algorithm with multiple 
subdomains for the model problem (1.1), and then analyze its convergence for the one 
dimensional wave equation. In Section 3 we present convergence result of the NNWR 
for multiple subdomains in 2D. Our convergence analysis shows that both the NNWR 
algorithm converge in a finite no of steps for finite time intervals, T < oo. Finally we 
present numerical results in Section 4, which illustrate our analysis. 

2. NNWR for multiple subdomains. In this section we define the Neumann- 
Neumann Waveform Relaxation (NNWR) method with many subdomains for the 
model problem (l.I) on the space-time domain D x (0,T) with Dirichlet data given 
on 5D. This can be treated as a generalization of the NNWR algorithm for two 
subdomains, for which see [12]. The method starts with a non-overlapping spatial 
domain decomposition, and the iteration involves subdomain solves in space time 
with corresponding interface condition, followed by a correction step. 

2.1. NNWR algorithm. Suppose the spatial domain D is partitioned into N 
non-overlapping subdomains {Di, 1 < i < N}, as illustrated in the left panel of 
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Figure 2.1. For i = 1,N set Ti := Dili \ Oil, Ai := {j S {1,..., TV} : Fi fl 
Fj has nonzero measure} and Fy := dQi fl dQj, so that the interface of Hi can be 
rewritten as F^ = UjeAi - We denote by Ui the restriction of the solution u of (1.1) 
to Hi and by riy the unit outward normal for Hi on the interface Fy . The NNWR 
algorithm for the model problem (1.1) starts with an initial guess w°^ (a;,t) along the 
interfaces Fi^ x (0, T), j S A^, * = 1,..., iV, and then performs the following two-step 
iteration: one first solves Dirichlet subproblems on each Hi in parallel, 


c^{x)Aui 

= /, 

in Hi, 

uf{x,0) 

= Uo{x), 

in H^, 

dtu’^{x,0) 

= Vo{x), 

in H^, 


= a, 

on dHi\Ti, 


= > 4 ''. 

on Tij,j e Ai 


One then solves Neumann subproblems on all subdomains, 


dtt^Pi - c^{x)A(p'y 
dt^p'^ ix,0) 


= 


with the updating condition 


0 , 

0 , 

0 , 

0 , 


> 


in Hi, 
in Oi, 
in Oi, 

on dHi \ Fj, 
on Fy ,j e Ai. 


wtj{x,t) = w'}- ^{x,t) - 9 {(fii |ryx(o,T) + |ryx(o,T)), 


( 2 . 1 ) 


( 2 . 2 ) 


(2.3) 


where 9 G (0,1] is a relaxation parameter. 

2.2. Convergence analysis for ID. We prove our convergence result for the 
one dimensional wave equation with constant speed, c{x) = c on the domain H := 
(0,L) with boundary conditions u{0,t) = go{t) and u{L,t) = gL{t), which in turn 
become zeros as we consider the error equations, f{x,t) = 0,go{t) = gL{t) = 0 = 
uo(x) = vo(x). We decompose H into non-overlapping subdomains := {xi-i,Xi), 
i = as shown in the right panel of Figure 2.1, and define the subdomain 

length hi := Xi — Xi-i, and hmin := mini<i<Ar Our initial guess is denoted by 
{^^(t)} on the interfaces {x = a;^} x (0, T), and for sake of consistency we denote 
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u>Q(i) = w^{t) = 0 for all k. We then obtain 


dttu'l - c^dxxUi = 0 , 


in rii, 




in Hi, dtt^Pi - c^dxx^i = 0 , 

w,^(x, 0 ) = 0 , inffi, V5*^(a:,0) = 0 , 

dtu^{x,0) = 0 , in rii, dtiPi{x,0) = 0 , 

u^{xi-i,t) = wfzlit), = idxUi_-^-dxu'^){xi-i,t), 


in rii, 




\t), 


dx(Piixi,t) = {dxu'l; - ^xu'|;_^_^){xi,t), 


(2.4) 

except for the first and last subdomains, where the Neumann conditions in the Neu¬ 
mann step are replaced by homogeneous Dirichlet conditions along the physical bound¬ 
aries. The updated interface values for the next step will be 


(t) = Wi \t)-0 {ip’y (xi, t) + {xi ,t)) . 


(2.5) 


We start by applying the Laplace transform to the homogeneous Dirichlet subproblems 
in (2.4), and obtain 

s'^Ui - C^Ui^xx = 0, Ui{Xi-i,s) = Wi-i{s), Ui{Xi,s) = Wi{s), 

for i = 2, ...,7V — 1. Set 7 i = cosh(/lis/c) ,ai = smh{his/c). These subdomain prob¬ 
lems have the solutions 

Ui{x, s) = — {wi{s) sinh {{x — Xi-i)s/c) + Wi-i(s) sinh {{xi — x)s/c)). 

Ui 

Next we apply the Laplace transform to the Neumann subproblems in (2.4) for sub- 
domains not touching the physical boundary, and obtain 

(pi{x, s) = Ci(s) cosh ((a; — Xi-i)s/c) -\- Di{s) cosh {{xi — x)s/c ), 

where 


C, = 

D, = 


'. / 7z 

l^+l\ 

Wi_i W,+ i 

Wi[ -h 



, \^i 


O'i <Ji+l 


, 7*-A 

Wi-2 Wi 

Wi-l ( — 

+ 


\<J^ 

O-j-l/ 

f-'i —1 


We therefore obtain for * = 2,..., — 2, at iteration k 

wHs) = wf"^(s) - 0 {^'l{xi,s)+lfi_^_J^{Xi,s)) 

= ^{s) — 9 {Ciji + Di + Ci+i + Z?i-|_i7i_|_i). 

Using the identity 7 ^^ — 1 = af and simplifying, we get 


w’^ = - , 


.fc-i 


27»7i-H 


^fe-l 


\ ^ ( 

1^+2 

7z\ 


J cr*+i V 

O'i+2 

UiJ 


7*-i l^+i 


^f+2 

^^2 




Cf i-\.2 Cri<7i-1 


■ ( 2 . 6 ) 
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For i = 1 and i = N, the Neumann conditions on the physical boundary are replaced 
by homogeneous Dirichlet conditions = 0 and (pN{L,t) = 0, t > 0. For these 

two subdomains, we obtain as solution after a Laplace transform 


(pi{x,s) = — 1 ( — + 

(Pn{x, s) = — ( WN-i 
JN 



\-k 

1 sinh ((a: 

cr2j 

02 ) 

'iN-l 1 7iv\ 

WN-2 

^crAT-i aN J 

<xn-i 


sinh {{xN — x)s/c). 


and thus the recurrence relations on the first interface is 
w'l = - 9 { w\-^ f 2 + 


7172 

1 cri72^ 

1 1 

f 73 

CTl\ 


(Ticr2 

7lO'2y' 

’ ' ^2 

VCTs 

7iy 

0'2U3 / 


and on the last interface, we obtain 


(2.7) 


Vn-I = w 


Nfe-1 

N-l 


- 9 


~fc-l 

^N-1 


iN-llN CTnIN-I 
O'N-lO'N InCTN-I 


+ - 


w 


N-2 


aN-1 


/ 7Af-2 

o-n\ 


\<JN-2 

IN J 

O'iV-lO'Ar-2 / 


( 2 . 8 ) 


We have the following convergence result for NNWR in ID: 

Theorem 2.1 (Convergence of NNWR for multiple subdomains). Let 9 = 1/4. 
Then the NNWR algorithm (2.4)-(2.5) converges in k+1 iterations, if the time window 
length T satisfies T/k < 2hminlc, c being the wave speed. 

Proof. With 0 = 1/4 the updating condition (2.6) becomes 


'(s) = -^ {kiw- ^(s) + ki+iwkkk + 

\k-l/ 


ti,i+2'ihj_|_2 (s) ii^i-2W^_2 (s)) , ( 2 . 9 ) 

where we defined k.i ■= {nn+i - <XiCTi+i), ii,i+i := := 

I _ 

-. Td 7_LO .- - 

iCTi+l 


■= FinW’ •= Similarly, we obtain for 

(2.7) 


vlis) = -i {kiwtkk + kawtkk-ks^tkk), 


( 2 . 10 ) 


where we defined ti i := (^ + ^-2), ii 2 = d_ ('t:! _ and £i 3 = 

V^ 7 i<T 2 0-1 (T 2 /’ 0-2 \^(T 3 71 / '’■ 2 Cr 3 

From (2.8), we obtain 


1 


w 


N 


_l(s) = —- (iN-l,N-lW%kk) + tN-l,N- 2 W%kk) ~ iN-l,N- 3 W%kk)) j 


where we dehned tAr_iAr_i = 


( div-l 
7N-1 


-l7iV-2 _|_ 7Af-l7iV-2 


-l(TN-2 (TN-lCrN-2 


— kj , i]S[-l,N-3 = 


( 2 . 11 ) 

1 


(TN-2(^N-3 


and tAr_i,Ar _2 = “ F^)' ^*4+2 ^ ii+ 2 ,^,ki+l = -U+l,i. So 


<7N-2 \<TN-3 in- 

by induction on (2.9)-(2.10) we can write 


w: 


i{s)— ^ ^ ( 4 ) Pi+j{ii+j,i+j-'iki+j,i+j — lT--M,i^---fi+j,i+j+lki+i,i+j+ 2 )w^_^j (s), 

j=-2n 


( 2 . 12 ) 
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and 


2n 

w1{s) = Y^ 

j=0 

where the coefficients are homogeneous polynomials of degree n. A similar ex¬ 
pression holds for w^_^{s). Now expanding hyperbolic functions into infinite binomial 
series, we obtain 


Pi+j [ii 


> ^i+j,2+j j ^i-i-i,3+i) 


w 


' k—Vi 

1+i 




(2.13) 


. = 2 cosh ((/l^ - h^+l)s/c) ^ ^ ^^- 2 his/c ^- 2 hi+,s/c'j 


sinh(his/c) sinh(ft,i_|_is/c) 




,—2himsjc 


oo oo 


^^ ^-2hi+ins{c _|_ EE^ 


m—l 

-2{mhi-\-nhi+i) s / c 


- _ _ sinh((fe; - hi+2)s/c) _ 

sinh{hiS/c) sinh(ft,i_|_is/c) sinh(ft,i_|_2s/c) 

OO OO OO 

E 2 "«'*t +2 X—^ X—^ 

e ^ + EE 


= 4 


OO oo 

2lsh^ ^^ 2msh^j^l 

e “ -f > e = 


E 2lsh^ V—> 

e “ + 


OO OO - 

+2 ^ ^ I 2{mhi+nhij^^)s 2{mh^j^i +nfe^^2) 

e <= -b e 


m—ln—1 
oo oo oo 

+EEE^ 

m—ln—1 


m—l 


2(Jh..j+Tnh..j_|_j^+n/ij_|_2)s 


^ — e ^ 


tiA+2 — 


1 


sinh(hi+is/c) sinh(ft,i+ 2 s/c) 


= 4e“(^*+i+^‘+2)'*/° 


i + E< 




_j_^^g-2nsfti+2/c 


OO oo 

EE‘ 


—2(m/^^^-l^-n/l.^+2)s/c 


^1.1 — 


2 cosh ((2/ii — h 2 )s/c) 
sinh(2/iis/c) sinh(ft, 2 s/c) 


_ ^ ^g-4his/c _|_ g-2/i2s/c^ 


1 + E' 


, — 4ms/^l / C 


OO OO 


_Ee-2"«'*2/c EE‘ 


— 2(2m/^l^-^^/^2)s/c 


n—1 


m—ln—1 


tl,2 = 


cosh ((hi — h^)s/c) 


cosh(his/c) sinh(h 2 s/c) sinh(h 3 s/c) 

OO OO OO 

n=l 


= 4 


l + ^(-l)'e-^ + ^< 


1=1 


2{mh-i+nh2)s , ^ _ 2 (mh,i+n/i3)s 2{mh2+rih2,)s 

-ij^e s- -b (-1) e s--b e ^- 


m—ln—1 
oo oo oo 


-EEE(-i) 

m—ln—1 


I ^ — 2{lhi-\-mh2-\-nh^)s/c 




+ h 2 )s/c _|_ ^-{h 2 + 2 h 3 )s/c 


)■ 
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The argument also holds similarly for the terms ^ 1 , 3 : tN-i,N- 2 , 

tN'-i,N- 3 - Now using these expressions we can write (2.12)-(2.13) as 


= (-1) 

and 


2k 


w 


lis) = i-lf 


(g-4/., 


j=-2/c 


2k 


, (2.14) 


—Ahiks/c _|_ — 2 h 2 ksj' 


=) w? (s) + ^ (s) (s) 


j=o 


(2.15) 


where (/^^^(s) and rj'^^ (s) are linear combinations of terms of the form e with 
m > 2khi/c for some I G {1,2,.. .N}. A similar expression holds for w%_i{s). We 
now recall the shifting property of Laplace transform 


where H{t) := 


£-1 I 

d, t > 0, 


^0, t<0. 

transform (2.14)-(2.15) and obtain 


{e-‘^^hs)}=H{t-a)f{t-a), (2.16) 

is the Heaviside step function. We use (2.16) to back 


= (t) = (-1)^ 


wnt- 


2khj 


H t- 


2khi 


+ wnt- 


2kh 


i+l 


H t- 


2khi. 


+1 


other terms], 




w. 


t - 


Akhi 


, ^kh\ » r| 

H[t -^ 1 + wl 


+ other terms] 


and a similar expression for So for T < 2kh^{^/c, we get wf(t) = 0 for all 

i, and the conclusion follows. □ 

Remark 2.2. The shifting property of Laplace transform (2.16) is the reason 
behind the finite step convergence of the DNWR for a particular value of the parameter 
9. The right hand side of (2.16) becomes identically zero for t < a, so that for 
sufficiently small time window length T (e.g., T < a) the error becomes zero and 
leads to convergence in the next iteration. In Figure 2.2 we plot C~^ 
fft) = sin(t) on the left, and show the effect of time-shifting on the right. 

3. Analysis of NNWR algorithm in 2D. We now formulate and analyze the 
NNWR algorithm for the two-dimensional wave equation 

dttu - c^Au = f{x, y, t), {x, y) G Ll = {I, L) x (0, tt), tG (0, T] 


with initial condition u{x, y, 0) = uo{x, y), dtu{x, y, 0) = vq{x, y) and Dirichlet bound¬ 
ary conditions. To define the Neumann-Neumann algorithm, we decompose fl into 
strips of the form Lli = (xi-i,Xi) x (0,7r), I = xq < Xi <■■■ < Xn = L. We define 
the subdomain width hi := Xi — Xi-i, and /imin := mini<i<ArAlso we directly 
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Fig. 2.2: Example of time-shifting for the function fit) = sin(t): C ^ 
left, and |e“®/(s)| on the right 


consider the error equations with f(x,y,t) = 0,uo{x,y) = 0 = vo{x,y) and homoge¬ 
neous Dirichlet boundary conditions. Given initial guesses {wi {y, t)} along the 

interface {x = Xi), the NNWR algorithm, as a particular case of (2.1)-(2.2)-(2.3), is 
given by performing iteratively for A: = 1, 2,... and for i = 1,... ,N the Dirichlet and 
Neumann steps 


dttUi - = 0, in 

wf(a:,2/,0) = 0, 

dtu'y{x,y,0) = 0 , 

= w^zl{y,t), 

u^{xt,y,t) = w'l~^{y,t), 

= u^{x,Tr,t) = 0, 


dttfi - c^^(Pi = 0, in Di, 

V5f(a;,J/,0) = 0, 
dty}’l{x,y,Q) = 0 , 

+ dniUf){Xi-i,y,t), 

dni'fHxi,y,t) = {dn,_^uf_i + dn,u'l){xi,y,t), 


(3.1) 

except for the first and last subdomain, where in the Neumann step the Neumann 
conditions are replaced by homogeneous Dirichlet conditions along the physical bound¬ 
aries. The update conditions are defined as 


{y,t) = Wi ^{y,t)-e{(pf{xi,y,t) + {xi,y,t)) . 


We perform a Fourier transform along the y direction to reduce the original problem 
into a collection of one-dimensional problems. Using a Fourier sine series along the 
y-direction, we get 


(a;, y,t) = Y^ U^ix, n, t) sin(ny) 

n>l 


where 


2 

Uf{x,n,t) = — / u^{x,r],t)sm{nri)dr]. 

TT Jo 


The NNWR algorithm (3.1) therefore becomes a sequence of ID problems for each n, 
’ -{x, n, t) — (? ^ (x, n, t) + c^n^Ui{x, n, t) = 0, 


dt^ 


dx"^ 


(3.2) 
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with the boundary conditions for U^{x, n, t). We now define 

x(a, /3, t) := |exp [—(3 ^/|, Re(s) > 0, (3.3) 

with s being the Laplace variable. Before presenting the main convergence theorem, 
we prove the following auxiliary result . 

Lemma 3.1. We have the identity: 

[o, 0 <t<p, 

where 5 is the dirac delta function and Ji is the Bessel function of first order given 
by 


1 r 

Ji{z) = — / cos {z simp — If) dip. 

Jo 


Proof. Using the change of variable r = \/s^ + we write 

g-/3r ^ g-/3. _ _ g-/3r)^ 

From the table [27, p. 245] we get 

£-i{e-^«}=5(t-/3). 

Also on page 263 of [27] we find 


£- 




0 < t < /3. 


(3.4) 


(3.5) 


Subtracting (3.5) from (3.4) we obtain the expected inverse Laplace transform. □ 
Now we are ready to prove the convergence result for NNWR in 2D: 

Theorem 3.2 (Convergence of NNWR in 2D). Let 6 = 1/4. For T > 0 fixed, 
the NNWR algorithm (3.1) converges in k + \ iterations, if the time window length T 
satisfies T/k < 2h,^ialc, c being the wave speed. 

Proof. We take Laplace transforms in t of (3.2) to get 

and now treat each n as in the one-dimensional analysis in the proof of Theorem 2.1, 
where the recurrence relations (2.9), (2.10) and (2.11) of the form 

w'f (s) = 51 Ay ^ (s) («) (3-6) 

3 

now become for each n = 1, 2,... 

Wi{n, s) = 55 ^\/-I- W^{n, s). 

3 


(3.7) 
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The equation (3.6) is of the form (2.14)-(2.15), that means A\j\s) are linear com¬ 
bination of terms of the form for g > 2khi/c for some I € {1,2,.. .N}. There¬ 
fore the coefficients (\/-I- c^n^) are sum of exponential functions of the form 

^ > 2khi/c. Hence we use the definition of x in (3-3) fo take the 
inverse Laplace transform of (3.7), and obtain 

= '^'^x{cn,gm,3A) * W°{n,t), 

j m 

with g-mj > 2khniin/c. So it is straightforward that for t < 2khmin/c, Wtin,t)=0 
for each n, since the function x is zero there by Lemma 3.1. Therefore the interface 
functions w^{y,t), given by w^{y,t) = ^lTj^(n, <) sin(ny) are also zero for all i G 

n>l 

{l,...,iV — 1}. Hence one more iteration produces the desired solution on the entire 
domain. □ 

4. Numerical Experiments. We perform numerical experiments to see the 
convergence behavior of the NNWR algorithm with multiple subdomains for the model 
problem 


dttu = c^{x)dxxu, X G (0,5), t > 0, 

m(x,0) = 0, Ut{x,0) = 0, 0 < a; < 5, (4.1) 

u(0, t) = u(5, t) = t > 0, 

which is discretized using centered finite differences in both space and time on a grid 
with Ax = At = 2x10“^. We consider a decomposition of (0,5) into five unequal 
subdomains, whose widths hi are 0.6, 0.6, 0.5, 2.3,1 respectively, and take the initial 
guesses w^{t) = G (0,7^] for 1 < j < 4. Note that in some of the experiments 
below, the coefficient c(a;) will be spatially varying. This will allow us to study how 
spatially varying coefficients affect the performance of the NNWR, which have only 
been analyzed in the constant coefficient case. For the first experiment, we take the 
constant speed, c = 1. On the left panel of Figure 4.1, we show the convergence for 
different values of the parameter 0 for T = 8 , and on the right the results for the 
best parameter 9 = 1/4 for different time window length T. We observe two-step 
convergence for 0 = 1/4 for a sufficiently small time window T. Now we take the 
propagation speed, c{x) = (a: + l )/6 for the second experiment. On the left panel of 
Figure 4.2, we show the convergence for different values of the parameter 0 for T = 8 , 
and on the right the results for the best parameter 0 = 1/4 for different time window 
length T. 

Next we show an experiment for the NNWR algorithm in two dimension for the 
following model problem 


dttU- {dxxU + dyyu) = Q,u{x, j/,0) = xy{x- l){y - tt){5x - 2){Ax - 3),Ut{x,y,0) = 0, 

with homogeneous Dirichlet boundary conditions. We discretize the wave equation 
using the centered finite difference in both space and time (Leapfrog scheme) on a 
grid with Ax = 5x10“^, Ay = 16x10“^, At = 4x10“^. We decompose our domain 
H := (0,1) X (0, tt) into three non-overlapping subdomains Hi = (0,2/5) x (0, tt), H 2 = 
(2/5, 3/4) X (0, tt), H 3 = (3/4,1) x (0, tt). As initial guesses, we take w°(y, t) = t sin(y). 
In Figure 4.3 we plot the convergence curves for different values of the parameter 0 for 
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Fig. 4.1: Convergence of NNWR with various values of 0 for T = 8 on the left, and 
for various lengths T of the time window and 9 = 1/4 on the right 




Fig. 4.2: Convergence of NNWR for variable coefficients c{x) = (a: + l)/6 with various 
values of 6 for T = 8 on the left, and for various lengths T of the time window and 
9 = 1/4 on the right 


T = 2 on the left panel, and on the right the results for the best parameter 0 = 1/4 
for different time window length T. 

We compare in Figure 4.4 the performance of the NNWR and DNWR (see [12]) 
algorithms with the SWR algorithms with and without overlap. Here we consider the 
problem 


dttu-d:^xU = 0, a; G (-3,2),t > 0, 

u{x,0) = 0, ut{x,0) = xe~^, — 3 < X < 2, 

= —3e^t, u{2,t) = 2te“^, t > 0, 

and for the overlapping Schwarz variant we use an overlap of length 24Ax, where 
Ax = 1/50. For the DNWR, NNWR and non-overlapping SWR we consider a domain 
decomposition into two subdomains fli = (—3,0) and D 2 = (0,2). We observe that 
the DNWR and NNWR algorithms converge as fast as the Schwarz WR algorithms 
for smaller time windows T. Due to the local nature of the Dirichlet-to-Neumann 
operator in ID [8], SWR converges in a finite number of iterations just like DNWR 
and NNWR. In higher dimensions, however, non-overlapping SWR will no longer 
converge in a finite number of steps, but DNWR and NNWR will; see Figure 4.5. 
Table 4.1 gives a summary of the theoretical results from Section 2 and 3 and [12], to 
indicate the maximum number of iterations needed for the ID and 2D wave equation 
to converge to the exact solution. For the comparison result in 2D, we consider the 
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Fig. 4.3: Convergence of NNWR in 2D: curves for different values of 6* for T = 2 on 
the left, and for various time lengths T and 0 = 1/4 on the right 




Fig. 4.4: Comparison of DNWR, NNWR, and SWR for ID wave equation for T = 4 
on the left, and T = 10 on the right 


model problem 


dttu - {dxxU + dyyu) = 0, u{x, y,0) = 0 = Ut{x, y, 0), 

with Dirichlet boundary conditions u{0,y,t) = sm{y),u{l,y,t) = y{y — Tr)t^ and 
u{x,0,t) = 0 = u(x, 7r,t). We decompose our domain D := (0,1) x (0,7r) for the two 
subdomains experiment into Di = (0,3/5) x (0,7r) and D 2 = (3/5,1) x (0,7r), and for 
the three subdomains experiment into Di = (0,2/5) x (0,7r), Q 2 = (2/5,3/4) x (0,7r), 
D 3 = (3/4,1) X (0,7r). We take a random initial guess to start the iteration, and for 
the overlapping SWR we use an overlap of length 2Ax in all the experiments. We 
implement first order methods with one parameter in optimized SWR iterations; for 
more details see [a]. On the left panel of Figure 4.5 we plot the comparison curves for 
two subdomains, and the same for three subdomains on the right. 

Now we show a numerical experiment for the NNWR algorithm with different 
time grids for different subdomains and discontinuous wave speed across interfaces. 
We consider the model problem 

dttu - c^dxxU = 0, u{x, 0) = 0 = utix, 0), 

with Dirichlet boundary conditions u{0,t) = Suppose the spatial 

domain D := (0,6) is decomposed into three equal subdomains Vti,i = 1,2,3, and 
the random initial guesses are used to start the NNWR iteration. For the spatial 
discretization, we take a uniform mesh with size Ax = 1x10“^, and for the time 
discretization, we use non-uniform time grids Ati,i = 1,2,3, as given in Table 4.2. 
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Table 4.1: Comparison of steps needed for convergence for ID wave equation. 


Methods 

2 subdomains, ID 

Many subdomains, ID 

Many subdomains, 2D 

DNWR 

T < 2fchmin/c 

T < khmin/c 

T <C ^^min/c 

NNWR 

T < Akh^irilc 

T < 2fchmin/c 

T < 2fchmin/c 


o 

a5 


10 “ 


10 -^ 


10 " 



-DNWR for 0=1/2 


NNWR for 0=1/4 


■■■Classical SWR 


.OSWR no overlap 


.OSWR 01 theor. 

\ \ \ 

■•OSWR 01 numer. 


3579 11 13 15 17 19 


iteration 



-DNWR for 0=1/2 


NNWR for 0=1/4 

\ ***♦- 

■■■Classical SWR 


.OSWR no overlap 

\ 

.OSWR 01 theor. 


^OSWR 01 numer. 


3 5 7 . 9 . 11 13 15 17 


iteration 


Fig. 4.5: Comparison of DNWR, NNWR, and SWR for T = 2 in 2D for two subdo¬ 
mains on the left, and three subdomains on the right 


For the non-uniform mesh grid, boundary data is transmitted from one subdomain 
to a neighboring subdomain by introducing a suitable time projection. For two di¬ 
mensional problems, the interface is one dimensional. Using ideas of merge sort one 
can compute the projection with linear cost, see [9] and the references therein. Even 
for higher dimensional interfaces, such an algorithm with linear complexity is still 
possible, see [10]. In Figure 4.6 we show the non-uniform time steps for different 
subdomains. Figure 4.7 gives the convergence behavior of the NNWR algorithm for 
T = 2 with the same non-uniform time grids as in Table 4.2. 

Finally we raise the issue of scalability of the NNWR algorithm by giving some 
numerical examples for the wave equation. From Theorem 2.1, one can say that 
as long as h/T is constant, we expect identical convergence behavior of the NNWR 
algorithm. We plot in Figure 4.8 the convergence curves by doubling the number of 
subdomains and making the time window length half. One can therefore conclude 
that the NNWR algorithm is weakly scalable for the wave equation. 

5. Conclusions . We defined the NNWR algorithm for multiple subdomains 
and for general hyperbolic problems, and analyzed their convergence properties for 
the second order wave equation in ID. We showed using numerical experiments that 
for a particular choice of the relaxation parameter, more specifically for 9 = 1/4, 
convergence can be achieved in a finite number of steps. In fact, this algorithm can 
be used as a two-step method, choosing the time window lengh T small enough. 
We have also extended the NNWR algorithm for the second order wave equation 
in 2D, and analyzed its convergence properties. We have also shown using numerical 
experiments that among the DNWR (see [11]) and NNWR methods, NNWR converges 
faster. But in comparison to DNWR, the NNWR has to solve twice the number of 
subproblems (once for Dirichlet subproblems, and once for Neumann subproblems) 
on each subdomain at each iteration. Therefore the computational cost is almost 
double for the NNWR than for the DNWR algorithm at each step. However, we get 
better convergence behavior with the NNWR in terms of iteration numbers. Finally 
we presented a comparison of performences between the DNWR, NNWR and Schwarz 
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Table 4.2: Propagation speed and time steps for different subdomains. 



Sll 

il.2 

fls 

wave speed c 

1/4 

2 

1/2 

time grids Ati 

13 X 10"^ 

39 X 10-^ 

1 X 10-^ 



Fig. 4.6: Subdomains with non-uniform time steps 
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Fig. 4.7: Solution using NNWR method for the wave equation with non-uniform time 
steps for 0 = 1/4 for T = 2: Solution after 1st iteration on the left, and solution after 
2nd iteration on the right 



Fig. 4.8: Graphs for comparing scalability of NNWR method for the wave equation 
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WR methods, and showed that the DNWR and NNWR converge faster than optimized 
SWR at least for higher dimensions. 
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